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We show that the two recently proposed methods to compute Renyi entanglement entropies in 
the realm of determinant quantum Monte Carlo methods for fermions are in principle equivalent, 
but differ in sampling strategies. The analogy allows to formulate a numerically stable calculation 
of the entanglement spectrum at strong coupling. We demonstrate the approach by studying static 
and dynamical properties of the entanglement hamiltonian across the interaction driven quantum 
phase transition between a topological insulator and quantum antiferromagnet in the Kane-Mele 
Hubbard model. The formulation is not limited to fermion systems and can readily be adapted to 
world-line based simulations of bosonic systems. 
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I. INTRODUCTION 

Consider a bipartition of a Hilbert space of a many 
body system in a state described by a density matrix p. 
Tracing over the degrees of freedom of one partition de¬ 
fines a reduced density matrix. Its entropy provides a 
measure of the entanglement between the two partitions 
[1]. At zero temperature one generically expects the en¬ 
tanglement entropy to follow an area law [2] . Corrections 
to this law have the potential of revealing fundamental 
properties such as topological order |3H5] or the central 
charge for one-dimensional systems [6]. The logarithm 
of the reduced density matrix defines an entanglement 
Hamiltonian [7] , the study of which has spurred substan¬ 
tial research IBHH]. The notion that it contains funda¬ 
mental and universal information has emerged and has 
been critically discussed m- The aim of this article is 
to develop tools to study the properties on the entangle¬ 
ment Hamiltonian in the realm of quantum Monte Carlo 
(QMC) simulations for fermions. 

For fermonic systems the calculation of the Renyi en¬ 
tanglement entropy has followed two different routes. 
One method builds on a replica idea with sampling based 
on a swap move [mus]. This approach was initially pro¬ 
posed for spin systems mms] at zero and finite tem¬ 
peratures and then generalized to fermions in the realm 
of determinant m and continuous time m quantum 
Monte Carlo (QMC) methods. We will refer to this al¬ 
gorithm as the swap algorithm. The other approach put 
forward in Ref. m utilizes the fact that in auxiliary field 
algorithms m - which express the interacting system in 
terms of a sum of non-interacting problems, the density 
matrix can be formally written as a sum over gaussian 
operators m- We will refer to this algorithm as the 
gaussian approach. It is in principle simple to implement 
and allows for generalizations to compute entanglement 
spectra [22] . As pointed out in [T71 [22] it suffers from an 
exponential growth of fluctuations in the strong coupling 
limit and when the subsystem size is large. 

We will show that within the auxiliary field approach 


both methods are equivalent, and merely correspond to 
different ways of carrying out the sampling. Since the 
swap algorithm is more stable than the gaussian one, the 
equivalence of the two methods shows how to stabilize 
the gaussian algorithm. As a consequence we are able 
to formulate a stable QMC algorithm allowing a detailed 
study of the entanglement Hamiltonian for fermion sys¬ 
tems at strong coupling. 

Here we will demonstrate the validity of the approach 
by studying a previously not accessible parameter region 
of the Kane-Mele Hubbard model [23U28] . In particu¬ 
lar, we will concentrate on the correlation induced phase 
transition from a topological insulator to a quantum an¬ 
tiferromagnetic from the perspective of the entanglement 
spectrum both in the single particle and particle-hole sec¬ 
tors. 

The article is organized as follows. In the next section 
we will show the equivalence of the swap and gaussian 
algorithms. Section [TTI| will use this equivalence to refor¬ 
mulate the proposed evaluation of the entanglement spec¬ 
tra of Ref. [22] in a numerically stable manner. Before 
concluding in Sec. [V| we test our approach by studying 
the correlation driven quantum phase transition in the 
Kane-Mele Hubbard model from the perspective of the 
entanglement spectrum. 


II. EQUIVALENCE SWAP AND GAUSSIAN 
ALGORITHMS FOR THE RENYI ENTROPY. 

Here we will start with the swap algorithm formula¬ 
tion of the Renyi entropy and derive the gaussian 
algorithm of Ref. m- We consider a real space parti¬ 
tioning of the Hilbert space, PL = PLa^PLb- To compute 
the n^^ Renyi entropy, 

5n = -^lnTr«>l, 


( 1 ) 



2 



FIG. 1. Schematic view of the imaginary time propagation 
of the Hamiltonian H(r) for the case n = 4. The Hamiltonian 
vanishes in the non-shaded regions. 

with Pa = Tv'j-i^p and p the density matrix, we consider 
the replicated Hilbert space: 

ntot = nA® ® nf ■ ■ ■ 'P'p. (2) 

At n = 1, 1-Ltot reduces to the original Hilbert space, 
1-L = 1-La and the Hamiltonian we will consider 

reads 

H = (3) 


Here 



(5) 
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A schematic representation of this time evolution is given 
in Fig. With this construction, one will show that: 

( 6 ) 

The above follows from writing the trace 

[O] = (7) 

where A and B run over a complete set of orthonormal 
states of Ha and Hb respectively and by noting that: 

{A,bA)^... = 

{A, B^^'> 1^1, SF) H • (8) 

i = l 

i^r 


In the swap algorithm one expands the imaginary time 
propagation from (3 to n/3 ( i.e. r G [0,n/3]) and defines 
a time dependent Hamiltonian in the Hilbert space Htot 
as: 

n 

H{r) = Y. 0 [r — (r — l)/3] 0 [r^ — r] (4) 


Here 0?=i such that corresponds to 

i^r 

the Hamiltonian in Hilbert space Ha G) H%^ . One can 
now explicitly compute the trace in Eq. (§ by inserting 
a complete set of states in Htot between each replica so 
as to obtain: 


= E E (Tli, An-i, ^ 


= E {Al\pA\An){An\pA\An-l)---{A2\pA\A^) = Z^TrH^ftX 


(9) 


Note that the reduced density matrix (A|p^|A') = 

;|r is independent on the 

choice of the replica. The partition function Z of the 
original Hamiltonian can be written as: 


E = Aru, 





( 10 ) 


d corresponds to the number of states per site ( d = 4 
for the spin-1/2 Hubbard model) and Nb the number of 
sites in the partition B such that counts the 

number of states in the Hilbert space H = (8)r=i • 

Hence, the factor compensates the over count¬ 


ing when computing the partition function of the original 
Hamiltonian by tracing over Htot- Note again that the 
partition function does not depend on the specific choice 
of the replica. 

Thus, 


= (11) 
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1 
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Tf’Wtot 

g-/3A(^) 

• • • 

g-/3A(i) 

d—n{n—l)NB 


We are now in the position to compute the Renyi en- 
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tropy with auxiliary field quantum Monte Carlo meth¬ 
ods. Here, we will use the finite temperature algorithm 
[29l [30]. For a given replica, we can make use of the 
Trotter decomposition so as to write 

^ ^ g-Arf(^)/2g-ArJ?«g-Arf(^)/2 ^ ^^^2) 

r=l 

( 12 ) 

and the Hubbard Stratonovitch transformation 

g-Arff<r> (13) 

For each imaginary time and replica, we have a vector 
of Hubbard Stratonovitch fields, It is important to 

(r) 

remember that, by construction, the dimension of s)- ^ 
is identical to that of a single simulation at n = 1. For 


the Hubbard model, corresponds to a vector of length 
Na^Nb of Ising spins and we have used a transformation 
where the Ising field couples to the local density [31]. 
and are single particle operators which one can 

write as: 

f{r) ^ (14) 


Here, c is a vector of fermionic annihilation operators 
running over all single particle states of the Hilbert space 
l-Ltot- The imaginary time propagation now reads: 


L-r 

g{r) T =1 


(15) 




(r) 

s(-) 


s(-) 


where is a short hand notation for 

With the above, we can compute the Renyi entropy as: 


SfiCi).. 


. . . 




fjin) 

. 


^-n(n-l)VB 


(16) 


TtbaP^ corresponds to the ratio of two partition func¬ 
tions, defined on the same configuration space. Note that 
the symmetries which ensure the absence of sign problem 
for the original Hamiltonian can be used to prove the ab¬ 
sence of sign problem for the numerator. For the Kane- 
Mele Hubbard model we refer the reader to [26l|27l|32] 
for a proof of the absence of sign problem at half-band 


filling. The ratio in Eq. can be computed with the 
swap algorithm described in [16]. This approach used 
to compute the Renyi entropies corresponds to the one 
adopted for bosonic systems and recently generalized to 
fermions m- To show the equivalence to the gaussian 
algorithm proposed in Ref. m and further developed in 
Ref. [22] to access entanglement spectrum we can rewrite 
Eq. (j^) as: 


,s(^) 


(17) 


where 




Tl'-Htot 

fjin) 

]. 





. 



(18) 


and 


((O)). 


(l)...s(n) 



. . . 




fjM 

. 


^-n(n-l)VB 


(19) 


The probability distribution P{s^^\ • • • is sampled 

by carrying out n-independent simulations of the original 
Hamiltonian. Our task is now to show that ((0))s(i)...s(n) 


reduces to Grover’s form m for the calculation of the 
Renyi entropy. 
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A. The n — 2 case. 

At n = 2 one can follow a pedestrian path and compute 
the ratio of the two fermionic determinants. We will 
sketch the calculation under the assumption that 
factorizes into spin-up and spin-down components such 
that we can only concentrate on the orbital degrees of 
freedom. Let Pa be a {Na + 2^^) x {Na + ‘^Nb) matrix 
with 


{Pa) 




5i^j Hi ^ A 
0 otherwise 


( 20 ) 


Here i and j run over all the single particle Wannier 
states of the Hilbert space Pa ^Pb ^Pb'^ and i e A 
states that Wanier state i belongs to Pa- Clearly Pa 
is a projector, and we will define similar quantities Pb 
and Pb'- Note that Pa^Pb and Pb' are projectors on 
orthogonal spaces such that for example = 0. For 

a given spin sector with d = 2 the integration over the 


((O)), 


(1).S(2) 


= det 


Pa (2G 


( 2 ) 

|( 2 ) 


PaG 


( 1 ) 

to 


Since Pa^Pb and Pb' are orthogonal projectors, the 
above determinant reduces to the determinant of the 


Na X Na matrix det 


{G^a 


1)(gW 




(r) (r) 

where corresponds to the Green function G^(^) re¬ 
stricted to Wannier states within Pa- The above is noth¬ 
ing but the equation put forward by Grover [19]. 


B. The general case. 


To show the equivalence for the Renyi entropy one 

(r) 

notes that acts non trivially in the Hibert space 

Pa^Pb^- Hence, 




(25) 


i=l 


The same calculation which leads to Eq. §) gives 


Prn, 

where 




= Tr„, 


(26) 


'pA{s^A = Pr^, 




(27) 


Using the relation 


Prn, 


U 


(r) 

(r) 




A) 

‘'sM 


(28) 


fermionic degrees of freedom gives [29[ [33| : 


((0))s(l),s(2) 


det 

t + G 

(2) /■/'(I) 

S(2)^S(1) 


det 

1 A. 

^ ^ ^S(2) 

det 

1 + 

2^ — 2Nb 


( 21 ) 


In the above equation we have defined 

L-r 


( 22 ) 


r=l 


Since the equal time Green function [33] in each replica 
reads. 


= 




1 -1 


(23) 


we can see, after some algebra, that 


^5(2) “ ^l(i) + + Pb + Pb' 


I 

one obtains: 

{{0))sA)...s(r^) = Tr^^ 
with 

Pa{s^A = 


(24) 


Pa{s^^^) ■ ■ ■ Pa{s^A (29) 


(1)1 


P'I’-ui’-) 

^ (P) 

u (L 



s{^) 



ruu 1 

u\\ 


(30) 


Pa{s^^^) is an operator acting in Pa- For a fixed Hub¬ 
bard Stratonovitch configuration, is a single particle 
propagator such that Wick’s theorem applies. As pointed 
out in [19] it has a Gaussian representation uniquely de- 

(r) 

fined by the Green function given at the end of the 
previous sub-section. In particular: 


A{sG^) = det ^1 - 


-d^ In 


^{r) 

1A 




(31) 


where a is a vector of fermionic annihilation operators 
running over all single particle states of the Hilbert space 
Pa- Taking the trace over Pa gives: 

((0))s(l)...s(n) = (32) 

det (l - det + R - ^ 

\ r=l ~ ^A 

which is nothing but the general result of Ref. m- 
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III. ENTANGLEMENT SPECTRA 


In Ref. [22] we proposed to compute the entanglement 
spectrum by considering the replica time displaced cor¬ 
relation function: 


5'o(^e) = (0^(ee)0)a 
Tr^A 




-o 


[p\] 


(33) 


for an operator O G Ha- Here te and n are integers with 
te < n. Within the gaussian approach m we can use 
the representation of the reduced density matrix, 


Pa = ^P{.s)pa{s), (34) 

S 


introduce n replicas and obtain: 


5g(rE) = 




( 1 )...«(») 


P{s 


( 1 ) 


("))Tr„^ U(s(’^^) ■ ■ PA{s^^^+^^)d^PA{s^^^'^) ■ ■ ■ PA{s^^^)d 


E»(i)...»w P(sW,--- [,3^(sW)---;3^(s(i))] 


(35) 


Sampling over n-independent simulations generates con¬ 
figurations distributed according to 
such that in principle one can compute numerator 
and denominator within a single simulation to pro¬ 
vide an estimate of the replica time displaced correla¬ 
tion function. This approach works at weak coupling 
but fails in the strong coupling limit due to fluctua¬ 
tions. Essentially, one is sampling the wrong distribu¬ 
tion, and re-weighting with the factor 


[pa{s‘' ^^) • • • which accounts for correla¬ 

tions between the replicas. Since one can show that the 
later quantity is positive it was proposed in [22] to sam¬ 
ple directly, 

so as to access the strong coupling regime. 

One can achieve this by using the above presented 
mapping between the gaussian and replica methods. In 
fact in the extended Hilbert space, one will see that: 


5'o(^e) 




Ue + 1) Af 77^^^ r) 

Es(i)...«(") Pr-Htot 

^5(1) 



(36) 


such that: 

{OirE)0)A= A«^'^---s^"^)((0(rE)0)),a)...,(.) 

(37) 

with 


finite temperature quantum Monte-Carlo methods. The 
above formulation is however not restricted to fermions. 
In fact, it carries over to bosonic systems amenable to 
stochastic simulations within, for example, the stochas¬ 
tic series expansion algorithm [34] . 


P(s(i)...s(«)) 


and 


Prutot 

fjin) 



Es(i)...«(") Pr-Htot 

Urb-- 



{{0{TE)0))s(i)...siu) = (39) 



o' 

P'^ntot 




The above corresponds to a standard calculation of an 
imaginary time displaced correlation function in the ex¬ 
tended Hilbert space at temperature n[3 albeit with an 
imaginary time dependent Hamiltonian. This quantity 
can readily be implemented in standard auxiliary field 


IV. RESULTS 

To illustrate the fact that we are able to access 
the strong coupling regime, we consider the interaction 
driven quantum phase transition in the Kane-Mele Hub¬ 
bard model. The model is defined on the Honeycomb 
lattice. Using the spinor notation c\ = it reads 

Hkmu = E 4 [% + i Y E (^1 ^» “ ^) • 

(40) 

The hopping matrix takes non-vanishing values, —t, 
between nearest neighbors of the honeycomb lattice, 
i — j = ±^ 1 , ±(^ 2 , (see Fig. and the intrinsic spin- 
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FIG. 2. Honeycomb lattice. For an L x L lattice, we con¬ 
sider periodic boundaries: Ci+^ai = Ci and Ci+La 2 = The 
real space partitioning breaks translation symmetry in the a 2 
direction but not along ai. 


Wa = 16,/3t = 4,n = 8 



FIG. 4. Entanglement spectral function as a function of the 
Hubbard U. For each /c-point the sum rule J dujAf i{k, u) — 1 
holds. In the plot, we have normalized the peak hight to unity. 


A = 0.2t, L = 12, Wa = 16, WB = 8,/3t = 4 



Te 


FIG. 3. Entanglement single particle Green function. Here 
we concentrate on the time reversal symmetric momentum 
k = TV and orbital m = 1 lying on the boundary of subsystem 

A. 


orbit term is given by 


^ij — ^ 


(i-r)x(r-i) 

|( 2 -r)x(r-i)| 

0 


if i, j are n.n.n. 
otherwise 


(41) 



q 


EIG. 5. x-component of the spin-spin correlation function 
taken on the edge of subsystem A corresponding to orbital 
index m = 1. 


time displaced Green function, 

G^,m'ik,rE) = (42) 

a 


where r is the intermediate site involved in the next near¬ 
est neighbor (n.n.n.) hopping process from site i to j. At 
A = 0.2t the model shows a zero temperature phase tran¬ 
sition between a quantum spin Hall state and a quantum 
antiferromagnetic at Uc/t = 5.71(2) [35]. The quantum 
phase transition is well understood and belongs to the 
3D XY universality class. Here, we show that we can de¬ 
tect this phase transition in the entanglement spectrum. 
In the absence of interactions the entanglement Hamil¬ 
tonian is adiabatically linked to the original one such 
that both have the same topological properties unni]. 
Thereby the entanglement Hamiltonian corresponding to 
a real space partitioning of the system should show edge 
states. 

Fig. |3] shows the single particle entanglement replica 


The real space cut we consider is translationally invariant 
in the ai lattice direction. Thereby, k = k • ai is a good 
quantum number which we can use to classify the data. 
The label m is an orbital index running across the width, 
Wa, of the cut. In Fig.j^we consider a 12 x 12 lattice with 
n = 8 replicas, Wa = 16, 20 and inverse temperatures 
/3t = 4,6,8. Note that Wa Wb = 21/ such that at 
n = S our largest simulations have 960 sites at an effective 
inverse temperature n^t = 64. All our simulations are 
carried out at a finite imaginary time step Art = 0.1. 
In Fig. [^ we concentrate on the time reversal symmetric 
momentum k = tv and orbital corresponding to the edge 
of the cut, m = 1. Since particle hole symmetry is present 
in the model, the Dirac cone is pinned at the fermi energy. 
Thereby a signature of the topological phase, is a non- 
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decaying single particle entanglement Green function as 
a function of the replica time te- At U/t = b we have 
considered various temperatures and values of Wa- As 
apparent, as a function of increasing Wa and thereby 
decreasing Wb^ Gfi{7r,TE) decays more quickly. This 
may be assigned to edge-edge correlations across the B 
subsystem. The phase transition is triggered by the onset 
of magnetic correlations which at T = 0 develop long 
range order beyond Uc thereby breaking time reversal 
symmetry. As a consequence, enhancing the temperature 
will reduced the magnetic correlation length, stabilize the 
topological state and show a less pronounced decay in 

To obtain a better overview of the data, we can define 
an entanglement spectral function by analytical continu¬ 
ation of the replica time data: 


1 r p-TEUJ 

Gm,m{k,TE) = - J duJYE^^^m,m{k,Uj)- ( 43 ) 


To carry out this step, we have used the stochastic Max¬ 
imum Entropy approach [asiisz]. Our results are plotted 
in Fig.[^ As apparent below Uc/t = 5.71(2) we observe a 
single Dirac cone and beyond the phase transition a gap 
in the entanglement spectrum opens. 

The gap in the entanglement spectral function stems 
from the onset of spin-spin correlations. The equal time 
spin-spin correlations of the entanglement Hamiltonian 
can be computed from 


{SlmSA,m')A ^ 


Tr 


PlSi 


OX 

q,m^—q,m' 


Tr[p- 


(44) 


with ^ Ei, corresponding 

to a partial Fourier transformation of the x-component of 
the spin-operator. As apparent from Fig. [^a sharp peak 
at g = TT emerges beyond the transition at Uc/t = 5.71(2) 


V. CONCLUSION 

In this article we have shown that the two methods put 
forward to compute the n^^ Renyi entropies in determi¬ 
nant QMC methods for fermions are in essence identical. 
Starting with the replica scheme proposed in [nms] and 
adapted to determinant m and continuous time m 
QMC, we can derive the free fermion or gaussian ap¬ 
proach put forward by Grover [19] . The two methods dif¬ 
fer in the sampling strategy. The gaussian approach sam¬ 
ples n independent replicas and correlations between the 
replicas are taken into account by re-weighting. The swap 
algorithm formulates the QMC in an extended Hilbert 
space thereby explicitly sampling correlations between 
replicas. The mapping between both methods shows 
how to formulate numerical simulations to access en¬ 
tanglement spectra at strong coupling by carrying out a 
standard simulation within the extended Hilbert space of 
subsystem A and n replicas of subsystem B albeit with 
a time dependent Hamiltonian. In contrast to our for¬ 
mer approach described in [22| the present formulation 
does not suffer from uncontrollable fluctuations in the 
strong coupling regime. We were able to study aspects 
of the entanglement spectrum in the correlation driven 
quantum phase transition between a topological insulator 
and quantum antiferromagnetic as realized in the Kane- 
Mele Hubbard model. The present formulation is numer¬ 
ically expensive since the total number of sites scales as 
Na + tiNb where n corresponds to the number of repli¬ 
cas and Na {Nb) the number of sites in subsystems A 
(B). The structure of the imaginary time evolution al¬ 
lows for many optimization strategies. Nevertheless, the 
overall computational effort scales as nP {Na tiNb)^^ 
Our approach to compute the entanglement spectrum 
is not specific to simulations of fermonic systems in the 
realm of determinant QMG methods. In fact it can be 
adapted to bosonic systems within, for example, the SSE 
[34] approach. Since these methods have a very favor¬ 
able scaling, nP {Na + uNb)^ introducing many replicas 
is not as expensive as for fermions. 
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